PrecipitationDaily.f90 Source File
Manage daily precipitation
Source Code
!! Manage daily precipitation !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.0 - 10th January 2019 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 10/Jan/2019 | Original code | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Module to manage daily precipitation ! MODULE PrecipitationDaily ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Parameters: float, short USE Loglib, ONLY : & !Imported routines: Catch USE ObservationalNetworks, ONLY : & !Imported types: ObservationalNetwork, & !Imported routines: ReadMetadata, ReadData, & CopyNetwork, & !Imported operators: ASSIGNMENT (=) USE Chronos, ONLY : & !Imported types: DateTime, & !Imported variables: timeString, & !Imported operations: OPERATOR (<), & OPERATOR (>), & OPERATOR (>=), & OPERATOR (<=), & OPERATOR (+), & OPERATOR (==), & ASSIGNMENT (=) USE GridLib, ONLY : & !Imported types: grid_real, grid_integer, & !Imported routines: NewGrid, GetDtGrid, & GridDestroy, ExportGrid, & SetLongName, SetStandardName, & SetUnits, SetCurrentTime, & !Imported parameters: NET_CDF, ESRI_ASCII, & ESRI_BINARY USE GridOperations, ONLY : & !Imported routines: GridByIni, GridResample, & CRSisEqual, GridConvert, & !Imported operations: ASSIGNMENT (=) USE IniLib, ONLY : & !Imported types: IniList, & !Imported routines: IniReadInt, IniReadReal, & KeyIsPresent, IniReadString, & SubSectionIsPresent USE Utilities, ONLY: & !Imported routines: GetUnit USE MeteoUtilities, ONLY: & !Imported routines: ReadField, ShiftMeteoWithLapse, & Interpolate, ShiftBackWithLapse, & Offset, Scale USE FileSys, ONLY: & !Imported routines: DirExists, FileExists, & FileDelete, CurrentDir, & GetOS, DirNew, & !Imported parameters: WIN32 USE GeoLib, ONLY: & !Imported variables: point1, point2, & !Imported routines: DecodeEPSG, Convert, & !Imported types: CRS, & !Imported operations: OPERATOR (==) USE StringManipulation, ONLY: & !Imported routines: ReplaceChar, ToString USE Units, ONLY: & !imported parameters: day IMPLICIT NONE ! Global (i.e. public) Declarations: INTEGER (KIND = short) :: dtPrecipitationDaily = 0 !!cumulation deltat TYPE(ObservationalNetwork) :: raingagesDaily !!raingage stations network TYPE(grid_real) :: precipitationRateDaily !![m/s] !Private declarations TYPE(DateTime), PRIVATE :: timeNewPrecipitation !!time when new precipitation data must be read TYPE(ObservationalNetwork), PRIVATE :: stationsRefElev !!stations network at reference elevation INTEGER (KIND = short), PRIVATE :: fileunit !! unit of file containing precipitation data INTEGER (KIND = short), PRIVATE :: interpolationMethod !!method to spatial interpolate site data INTEGER (KIND = short), PRIVATE :: interpolationMethod_assignment !!method to assign spatial interpolation ! !!1 = one method for the entire domain, 2 = a map with interpolation method codes TYPE (grid_integer), PRIVATE :: interpolationMethod_map INTEGER (KIND = short), PRIVATE :: interpolationMethod_vector (3) !!defines active interpolation methods TYPE (grid_real), PRIVATE :: interpolatedMap (3) !!1 = map for thiessen, 2 = map for idw, 3 = map for kriging REAL (KIND = float), PRIVATE :: cellsizeInterpolation !!spatial resolution of interpolated grid !TYPE(grid_real), PRIVATE :: gridTemp !temporary grid used to interpolate to ! !different spatial resolution !TYPE(grid_real), PRIVATE :: demTemp !temporary dem with the same resolution of precipitationTemp ! !used to apply lapse rate INTEGER (KIND = short), PRIVATE :: neighbors !!number of closest data to use for interpolation TYPE(grid_real), PRIVATE :: grid_devst !!standard deviation of kriging interpolation CHARACTER (LEN = 300), PRIVATE :: fileGrid !!file containing grids INTEGER (KIND = short), PRIVATE :: dtGrid !!dt of imported field INTEGER (KIND = short), PRIVATE :: elevationDrift !! 1 = use elevation to modify interpolated data !TYPE(grid_real), PRIVATE :: dem !!digital elevation model to be used to modify precipitation !TYPE(grid_real), PRIVATE :: dem_coarse !!digital elevation model at coarse resolution with ! !! CRS of hi-res dem !TYPE(grid_real), PRIVATE :: dem_coarse_temp !!temporary dem TYPE(grid_real), PRIVATE :: lapse_map !!lapse rate map !TYPE(grid_real), PRIVATE :: lapse_map_temp !!lapse rate map used for interpolation at coarse resolution REAL (KIND = float), PRIVATE :: lapse !!lapse rate [m/s/m] INTEGER (KIND = short), PRIVATE :: lapse_rate_assignment !!1 = scalar, 2 = map that may change with time INTEGER (KIND = short), PRIVATE :: i, j !!loop index REAL (KIND = float), PARAMETER, PRIVATE :: refElevation = 1000. !!reference elevation for applying lapse rate REAL (KIND = float), PRIVATE, PARAMETER :: MISSING_DEF_REAL = -9999.9 REAL (KIND = float), PRIVATE :: scale_factor !!scale factor to apply to precipitation map REAL (KIND = float), PRIVATE :: offset_value !!offset to apply to precipitation map REAL (KIND = float), PRIVATE :: valid_prcn !!when data from several time steps are read, !!this is the minimum percentage (0-1) of valid !!data that must be prresent to consider !!valid the aggregated value. !used by inverse distance weighting interpolation REAL (KIND = float), PRIVATE :: idw_power !! power to be used with IDW !used by kriging INTEGER (KIND = short), PRIVATE :: krige_var !! when set to 1 a grid of kriging variance is generated and exported if export option is active, default to 0 INTEGER (KIND = short), PRIVATE :: krige_anisotropy !! considers anisotropy, default = 0 excludes anisotropy INTEGER (KIND = short), PRIVATE :: krige_varmodel !! 1 = spherical, 2 = exponential, 3 = gaussian, 0 = automatic fitting. default to 2 INTEGER (KIND = short), PRIVATE :: krige_lags !! number of lags for semivariogram. if undefined or set to 0 default to 15 REAL (KIND = float) , PRIVATE :: krige_maxlag !! maximum distance to be considered for semivariogram assessment. If undefined or set to 0, it is computed automatically INTEGER (KIND = short), PRIVATE :: export !!activates grid exporting CHARACTER (LEN = 1000), PRIVATE :: export_path !!folder where to put exported grids CHARACTER (LEN = 1000), PRIVATE :: export_file !!name of exported precipitation grid CHARACTER (LEN = 1000), PRIVATE :: export_file_var !!name of exported precipitation variance grid TYPE (DateTime), PRIVATE :: export_start !! time and date to start exporting TYPE (DateTime), PRIVATE :: export_stop !! time and date to stop exporting INTEGER (KIND = short), PRIVATE :: export_dt !! time between two exportations TYPE (DateTime), PRIVATE :: timeNewExport !! time when new exporting must occur INTEGER (KIND = short), PRIVATE :: export_format !! 1 = esri_ascii, 2 = esri_binary, 3 = netcdf INTEGER (KIND = short), PRIVATE :: export_epsg !! coordinate reference system of exported grid TYPE (grid_real), PRIVATE :: exportedGrid, exportedGridVar TYPE (CRS), PRIVATE :: exportGridMapping LOGICAL, PRIVATE :: needConversion !Public routines PUBLIC :: PrecipitationDailyInit PUBLIC :: PrecipitationDailyRead !Local routines PRIVATE :: SetSpecificProperties !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize precipitation SUBROUTINE PrecipitationDailyInit & ! ( ini, mask, dtMeteo, tstart, dem_loaded ) IMPLICIT NONE TYPE (IniList), INTENT(IN) :: ini TYPE (grid_integer), INTENT(IN) :: mask !!defines interpolation extent INTEGER (KIND = short), INTENT(IN) :: dtMeteo !! deltat of meteo data reading TYPE (DateTime), INTENT(IN) :: tstart !!initial time LOGICAL , INTENT (in) :: dem_loaded !! true if dem has been loaded !local declarations CHARACTER (LEN = 1000) :: filename TYPE (DateTime) :: lapseTime TYPE (grid_real) :: meteoTemp !-------------------------end of declarations---------------------------------- !set initial time timeNewPrecipitation = tstart !initialize grid CALL NewGrid (precipitationRateDaily, mask, 0.) CALL NewGrid (grid_devst, mask, 0.) !set deltat dtPrecipitationDaily = day !check dtPrecipitation is multiple of dtMeteo IF (.NOT.(MOD(dtPrecipitationDaily,dtMeteo) == 0)) THEN CALL Catch ('error', 'PrecipitationDaily', & 'dt is not multiple of dtMeteo') END IF !set valid threshold IF (KeyIsPresent ('valid-threshold', ini, & section = 'precipitation-daily') ) THEN valid_prcn = IniReadReal ('valid-threshold', ini, & section = 'precipitation-daily') ELSE ! set to default = 1.0 CALL Catch ('info', 'PrecipitationDaily', & 'valid-threshold not defined, set to 1') valid_prcn = 1. END IF !set cell-size cellsizeInterpolation = mask % cellsize !interpolation-assignment method IF (KeyIsPresent('interpolation-assignment', ini, & section = 'precipitation-daily') ) THEN interpolationMethod_assignment = IniReadInt ('interpolation-assignment', ini, & section = 'precipitation-daily') ELSE CALL Catch ('error', 'PrecipitationDaily', & 'interpolation-assignment missing in meteo configuration file') END IF !set interpolation method IF (interpolationMethod_assignment == 1) THEN interpolationMethod = IniReadInt ('interpolation', ini, & section = 'precipitation-daily') CALL SetSpecificProperties (interpolationMethod, ini, mask) ELSE !read map interpolationMethod = -1 CALL GridByIni (ini, interpolationMethod_map, "precipitation-daily", & "interpolation") !scan for interpolation methods included interpolationMethod_vector = 0 DO i = 1, interpolationMethod_map % idim DO j = 1, interpolationMethod_map % jdim IF (interpolationMethod_map % mat(i,j) /= & interpolationMethod_map % nodata) THEN interpolationMethod_vector & (interpolationMethod_map % mat (i,j)) = & interpolationMethod_map % mat (i,j) END IF END DO END DO DO i = 1, 3 CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask) END DO END IF !scale factor and offset IF (KeyIsPresent('offset', ini, section = 'precipitation-daily')) THEN offset_value = IniReadReal ('offset', ini, section = 'precipitation-daily') ELSE offset_value = MISSING_DEF_REAL END IF IF (KeyIsPresent('scale-factor', ini, section = 'precipitation-daily')) THEN scale_factor = IniReadReal ('scale-factor', ini, section = 'precipitation-daily') ELSE scale_factor = MISSING_DEF_REAL END IF !set power_idw IF (KeyIsPresent('idw-power', ini, section = 'precipitation-daily')) THEN idw_power = IniReadReal ('idw-power', ini, section = 'precipitation-daily') ELSE !set default value idw_power = 2. END IF !file filename = IniReadString ('file', ini, section = 'precipitation-daily') IF (interpolationMethod_assignment == 1 .AND. & interpolationMethod == 0) THEN !data are stored in net-cdf file !store net-cdf filename fileGrid = filename IF ( KeyIsPresent ('variable', ini, section = 'precipitation-daily') ) THEN precipitationRateDaily % var_name = IniReadString ('variable', & ini, section = 'precipitation-daily') !read grid and store as temporary file CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), & fileFormat = NET_CDF, & variable = TRIM(precipitationRateDaily % var_name) ) ELSE IF (KeyIsPresent ('standard_name', ini, & section = 'precipitation-daily') ) THEN precipitationRateDaily % standard_name = IniReadString ('standard_name', & ini, section = 'precipitation-daily') ELSE CALL Catch ('error', 'PrecipitationDaily', & 'standard_name or variable missing in section precipitation-daily' ) END IF !set cellsize to zero. Optimal cellsize is automatically computed cellsizeInterpolation = 0. ELSE !open file containing local measurements fileunit = GetUnit () OPEN(fileunit,file = filename(1:LEN_TRIM(filename)),status='old') END IF !use elevation to modify precipitation IF ( KeyIsPresent ('elevation-drift', ini, section = 'precipitation-daily') ) THEN elevationDrift = IniReadInt ('elevation-drift', ini, section = 'precipitation-daily') ELSE elevationDrift = 0 !default, suppress drift END IF IF (elevationDrift == 1) THEN IF (interpolationMethod == 0 ) THEN CALL Catch ('error', 'PrecipitationDaily', & 'elevation drift cannot be applied when interpolation = 0 ') END IF !check if dem have been loaded by domain properties IF ( .NOT. dem_loaded) THEN CALL Catch ('error', 'PrecipitationDaily', & 'dem for elevation drift was not loaded ') END IF !set lapse rate assignment IF (KeyIsPresent ('lapse-rate-assignment', ini, & section = 'precipitation-daily') ) THEN lapse_rate_assignment = IniReadInt ('lapse-rate-assignment', ini, & section = 'precipitation-daily') ELSE !lapse-rate-assignment missing CALL Catch ('error', 'PrecipitationDaily', & 'lapse-rate-assignment missing in meteo configuration file') END IF !lapse rate map IF (lapse_rate_assignment == 2) THEN IF ( SubSectionIsPresent (subsection = 'lapse-rate-map', & section = 'precipitation-daily', iniDb = ini) ) THEN !check all keywords are defined IF (IniReadString ('format', ini, 'precipitation-daily',& 'lapse-rate-map') == 'net-cdf') THEN IF (KeyIsPresent ('time',ini, 'precipitation-daily', & 'lapse-rate-map')) THEN timeString = IniReadString ('time', ini, & 'precipitation-daily', 'lapse-rate-map') lapseTime = timeString IF (lapseTime > tstart) THEN CALL Catch ('error', 'PrecipitationDaily', & 'lapse-rate-map time greater than starting time' ) END IF ELSE CALL Catch ('error', 'PrecipitationDaily', & 'time keyword is missing in lapse-rate-map subsection' ) END IF END IF CALL GridByIni (ini, lapse_map, "precipitation-daily", "lapse-rate-map") !check coordinate reference system and spatial resolution IF (.NOT. CRSisEqual (mask = mask, grid = lapse_map, & checkCells= .TRUE.)) THEN CALL Catch ('error', 'PrecipitationDaily', & 'wrong lapse-rate-map spatial reference system ') END IF ELSE CALL Catch ('error', 'PrecipitationDaily', & 'lapse-rate-map missing in meteo configuration file') END IF ELSE !build lapse map from scalar value IF ( KeyIsPresent ('lapse-rate-scalar', ini, & section = 'precipitation-daily') ) THEN lapse = IniReadReal ('lapse-rate-scalar', ini, & section = 'precipitation-daily') CALL NewGrid (lapse_map, mask) lapse_map = lapse / 3600. !convert lapse rate from mm/h/m to mm/s/m ELSE CALL Catch ('error', 'PrecipitationDaily', & 'lapse-rate-scalar missing in meteo configuration file') END IF END IF !lapse_rate_assignment END IF !elevationDrift !grid exporting settings IF (KeyIsPresent ('export', ini, section = 'precipitation-daily') ) THEN export = IniReadInt ('export', ini, section = 'precipitation-daily') ELSE export = 0 END IF IF (export == 1) THEN !set export path name IF (KeyIsPresent ('export-path', ini, section = 'precipitation-daily') ) THEN export_path = IniReadString ('export-path', ini, section = 'precipitation-daily') !check if path ends with / or \ IF ( export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '\' .AND. & export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '/' ) THEN export_path (LEN_TRIM (export_path)+1:LEN_TRIM (export_path)+1) = '\' END IF IF (INDEX (export_path,'.') == 1) THEN !detected relative path, remove '.' export_path = export_path (2:LEN_TRIM(export_path)) !build absolute path export_path = TRIM(CurrentDir() ) // TRIM(export_path) END IF !check OS convention IF (GetOS () == WIN32) THEN export_path = ReplaceChar (export_path,'/','\') ELSE export_path = ReplaceChar (export_path,'\','/') END IF !check folder exists IF ( .NOT. DirExists (TRIM (export_path) ) ) THEN CALL Catch ('info', 'PrecipitationDaily', & 'creating directory: ', argument = TRIM(export_path)) CALL DirNew (export_path) END IF ELSE CALL Catch ('error', 'PrecipitationDaily', & 'missing export-path') END IF !starting time IF (KeyIsPresent ('export-start', ini, section = 'precipitation-daily') ) THEN timeString = IniReadString ('export-start', ini, 'precipitation-daily') export_start = timeString ELSE CALL Catch ('error', 'PrecipitationDaily', & 'missing export-start') END IF !initialize timeNewExport timeNewExport = export_start !ending time IF (KeyIsPresent ('export-stop', ini, section = 'precipitation-daily') ) THEN timeString = IniReadString ('export-stop', ini, 'precipitation-daily') export_stop = timeString ELSE CALL Catch ('error', 'PrecipitationDaily', & 'missing export-start') END IF !export dt IF (KeyIsPresent ('export-dt', ini, section = 'precipitation-daily') ) THEN export_dt = IniReadInt ('export-dt', ini, section = 'precipitation-daily') ELSE CALL Catch ('error', 'PrecipitationDaily', & 'missing export-dt') END IF !coordinate reference system IF (KeyIsPresent ('export-epsg', ini, section = 'precipitation-daily') ) THEN export_epsg = IniReadInt ('export-epsg', ini, section = 'precipitation-daily') exportGridMapping = DecodeEPSG (export_epsg) IF (exportGridMapping == precipitationRateDaily % grid_mapping) THEN needConversion = .FALSE. !initialize grid for exporting with CRS as precipitationRate CALL NewGrid (exportedGrid, precipitationRateDaily) CALL NewGrid (exportedGridVar, precipitationRateDaily) ELSE needConversion = .TRUE. !initialize grid for converting with required CRS exportedGrid % grid_mapping = DecodeEPSG (export_epsg) exportedGridVar % grid_mapping = DecodeEPSG (export_epsg) END IF ELSE CALL Catch ('info', 'PrecipitationDaily', & 'export-epsg not defined, use default') needConversion = .FALSE. !initialize grid for exporting with CRS as precipitationRate CALL NewGrid (exportedGrid, precipitationRateDaily) CALL NewGrid (exportedGridVar, precipitationRateDaily) END IF !export file format IF (KeyIsPresent ('export-format', ini, section = 'precipitation-daily') ) THEN export_format = IniReadInt ('export-format', ini, section = 'precipitation-daily') ELSE CALL Catch ('error', 'PrecipitationDaily', & 'missing export-format') END IF IF (export_format == 3) THEN !precipitation map CALL SetLongName ( 'precipitation_amount', exportedGrid) CALL SetStandardName ( 'precipitation_amount', exportedGrid) CALL SetUnits ('kg m-2', exportedGrid) !this unit is for exporting grid, it is converted internally to m/s !if file exists, remove it export_file = TRIM(export_path) // 'precipitation_daily.nc' IF ( FileExists (export_file) ) THEN CALL FileDelete (export_file) END IF !variance of kriging CALL SetLongName ( 'precipitation_amount_variance', exportedGridVar) CALL SetStandardName ( 'precipitation_amount_variance', exportedGridVar) CALL SetUnits ('kg m-2', exportedGridVar) !this unit is for exporting grid !if file exists, remove it export_file_var = TRIM(export_path) // 'precipitation_daily_variance.nc' IF ( FileExists (export_file_var) ) THEN CALL FileDelete (export_file_var) END IF END IF END IF !complete initialization IF (interpolationMethod == 0) THEN !Get the dt of imported field. Assume dt is regular dtGrid = GetDtGrid (filename = fileGrid, checkRegular = .TRUE.) !check dt is multiple of dtGrid IF (.NOT.(MOD(dtPrecipitationDaily,dtGrid) == 0)) THEN CALL Catch ('error', 'PrecipitationDaily', & 'dt precipitation is not multiple of dt of input grid') END IF ELSE !populate raingages metadata CALL ReadMetadata (fileunit, raingagesDaily) !check spatial reference system IF ( .NOT. raingagesDaily % mapping == mask % grid_mapping) THEN CALL Catch ('info', 'PrecipitationDaily', & 'converting coordinate of stations') !convert stations' coordinate point1 % system = raingagesDaily % mapping point2 % system = mask % grid_mapping raingagesDaily % mapping = mask % grid_mapping DO i = 1, raingagesDaily % countObs point1 % easting = raingagesDaily % obs (i) % xyz % easting point1 % northing = raingagesDaily % obs (i) % xyz % northing point1 % elevation = raingagesDaily % obs (i) % z CALL Convert (point1, point2) raingagesDaily % obs (i) % xyz % easting = point2 % easting raingagesDaily % obs (i) % xyz % northing = point2 % northing END DO END IF !set supplementary raingages network for elevation drift IF (elevationDrift == 1) THEN stationsRefElev = raingagesDaily DO i = 1, raingagesDaily % countObs stationsRefElev % obs (i) % xyz % elevation = refElevation END DO END IF END IF RETURN END SUBROUTINE PrecipitationDailyInit !============================================================================== !| Description: ! Read precipitation data SUBROUTINE PrecipitationDailyRead & ! ( time, dem) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !!current time TYPE (grid_real), INTENT(IN) :: dem !local declarations: TYPE (DateTime) :: time_toread !! time to start reading from CHARACTER (LEN = 300) :: filename CHARACTER (LEN = 300) :: varname INTEGER (KIND = short) :: i, j !-------------------------end of declarations---------------------------------- IF ( .NOT. (time < timeNewPrecipitation) ) THEN !time_toread = time + - (dtPrecipitationDaily - raingagesDaily % timeIncrement) time_toread = time timeString = time CALL Catch ('info', 'PrecipitationDaily', & 'read new precipitation data: ', & argument = timeString) !update lapse rate when required IF (elevationDrift == 1) THEN IF (time == lapse_map % next_time) THEN !update lapse_map varname = lapse_map %var_name filename = lapse_map % file_name !destroy current grid CALL GridDestroy (lapse_map) !read new grid in netcdf file CALL NewGrid (layer = lapse_map, fileName = TRIM(filename), & fileFormat = NET_CDF, & variable = TRIM(varname), & time = time_toread) END IF END IF SELECT CASE (interpolationMethod) CASE (0) !input grid CALL ReadField (fileGrid, time_toread, & dtPrecipitationDaily, dtGrid, & 'C', precipitationRateDaily, & varName = precipitationRateDaily % var_name) CASE DEFAULT !requires interpolation !read new precipitation data CALL ReadData (network = raingagesDaily, fileunit = fileunit, & time = time_toread, aggr_time = dtPrecipitationDaily, & aggr_type = 'sum', tresh = valid_prcn) IF (elevationDrift == 1) THEN !shift precipitation observation to reference elevation by applying lapse rate CALL ShiftMeteoWithLapse (raingagesDaily, lapse_map, refElevation, & stationsRefElev, dtPrecipitationDaily) !interpolate IF (interpolationMethod_assignment == 1) THEN !only one method is applied CALL Interpolate (network = stationsRefElev, & method = interpolationMethod, & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = precipitationRateDaily, & devst = grid_devst) ELSE !loop trough interpolation methods DO i = 1, 3 IF (interpolationMethod_vector(i) > 0) THEN CALL Interpolate (network = stationsRefElev, & method = interpolationMethod_vector(i), & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = interpolatedMap (i), & devst = grid_devst) END IF END DO !compose the final interpolated field DO i = 1, interpolationMethod_map % idim DO j = 1, interpolationMethod_map % jdim IF (interpolationMethod_map % mat(i,j) /= & interpolationMethod_map % nodata ) THEN precipitationRateDaily % mat (i,j) = & interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j) END IF END DO END DO END IF !Shift back interpolated field to terrain surface CALL ShiftBackWithLapse (precipitationRateDaily, dem, & lapse_map, refElevation, & dtPrecipitationDaily) ELSE IF (interpolationMethod_assignment == 1) THEN !only one method is applied CALL Interpolate (network = raingagesDaily, & method = interpolationMethod, & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = precipitationRateDaily, & devst = grid_devst) ELSE !loop trough interpolation methods DO i = 1, 3 IF (interpolationMethod_vector(i) > 0) THEN CALL Interpolate (network = raingagesDaily, & method = interpolationMethod_vector(i), & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = interpolatedMap (i), & devst = grid_devst) END IF END DO !compose the final interpolated field DO i = 1, interpolationMethod_map % idim DO j = 1, interpolationMethod_map % jdim IF (interpolationMethod_map % mat(i,j) /= & interpolationMethod_map % nodata ) THEN precipitationRateDaily % mat (i,j) = & interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j) END IF END DO END DO END IF !1 or more interpolation methods END IF !elevationDrift END SELECT !apply scale factor and offset, if defined IF (offset_value /= MISSING_DEF_REAL) THEN CALL Offset (precipitationRateDaily, offset_value) END IF IF (scale_factor /= MISSING_DEF_REAL) THEN CALL Scale (precipitationRateDaily, scale_factor) END IF !filter data < 0. DO i = 1, precipitationRateDaily % idim DO j = 1, precipitationRateDaily % jdim IF (precipitationRateDaily % mat(i,j) /= precipitationRateDaily % nodata) THEN IF ( precipitationRateDaily % mat(i,j) < 0.) precipitationRateDaily % mat(i,j) = 0. END IF END DO END DO !grid exporting IF(export > 0 ) THEN IF( time == timeNewExport .AND. & time >= export_start .AND. & time <= export_stop ) THEN timeString = time timeString = timeString(1:19) timeString(14:14) = '-' timeString(17:17) = '-' grid_devst % reference_time = precipitationRateDaily % reference_time IF (needConversion) THEN CALL GridConvert (precipitationRateDaily, exportedGrid) CALL GridConvert (grid_devst, exportedGridVar) ELSE exportedGrid = precipitationRateDaily exportedGridVar = grid_devst END IF SELECT CASE (export_format) CASE (1) !esri-ascii export_file = TRIM(export_path) // TRIM(timeString) // & '_precipitation.asc' CALL Catch ('info', 'PrecipitationDaily', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII) IF (krige_var == 1) THEN !export kriging variance export_file_var = TRIM(export_path) // TRIM(timeString) // & '_precipitation_variance.asc' CALL Catch ('info', 'PrecipitationDaily', & 'exporting variance grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII) END IF CASE (2) !esri-binary export_file = TRIM(export_path) // TRIM(timeString) // & '_precipitation' CALL Catch ('info', 'PrecipitationDaily', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY) IF (krige_var == 1) THEN !export kriging variance export_file_var = TRIM(export_path) // TRIM(timeString) // & '_precipitation_variance' CALL Catch ('info', 'PrecipitationDaily', & 'exporting variance grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY) END IF CASE (3) !net_cdf CALL SetCurrentTime (time, exportedGrid) CALL Catch ('info', 'PrecipitationDaily', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, 'precipitation_amount', 'append') IF (krige_var == 1) THEN !export kriging variance CALL SetCurrentTime (time, exportedGridVar) CALL Catch ('info', 'PrecipitationDaily', & 'exporting grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, 'precipitation_amount_variance', 'append') END IF END SELECT timeNewExport = timeNewExport + export_dt END IF ENDIF !conversion mm => m/s DO i = 1, precipitationRateDaily % idim DO j = 1, precipitationRateDaily % jdim IF (precipitationRateDaily % mat(i,j) /= precipitationRateDaily % nodata) THEN precipitationRateDaily % mat(i,j) = precipitationRateDaily % mat(i,j) / & dtPrecipitationDaily / 1000. END IF END DO END DO !time forward timeNewPrecipitation = timeNewPrecipitation + dtPrecipitationDaily END IF RETURN END SUBROUTINE PrecipitationDailyRead !============================================================================== !| Description: ! set properties and initialize variables for each interpolation method SUBROUTINE SetSpecificProperties & ! ( method, ini, mask ) IMPLICIT NONE !Arguments with intent(in): INTEGER (KIND = short), INTENT(in) :: method TYPE (IniList) , INTENT(in) :: ini TYPE (grid_integer) , INTENT(in) :: mask !----------------------------end of declarations------------------------------- !set data specific for each interpolation method SELECT CASE (method) CASE (1) !thiessen CALL NewGrid (interpolatedMap (1), mask, 0.) CASE (2) !inverse distance CALL NewGrid (interpolatedMap (2), mask, 0.) !nearest-points IF (KeyIsPresent ('nearest-points', ini, & section = 'precipitation-daily') ) THEN neighbors = IniReadInt ('nearest-points', ini, & section = 'precipitation-daily') ELSE !nearest-points missing CALL Catch ('error', 'PrecipitationDaily', & 'nearest-points missing in meteo configuration file') END IF CASE (3) !kriging CALL NewGrid (interpolatedMap (3), mask, 0.) !nearest-points IF (KeyIsPresent ('nearest-points', ini, & section = 'precipitation-daily') ) THEN neighbors = IniReadInt ('nearest-points', ini, & section = 'precipitation-daily') CALL Catch ('info', 'PrecipitationDaily', 'neighbors set to: ', argument = ToString (neighbors) ) ELSE !nearest-points missing CALL Catch ('error', 'PrecipitationDaily', & 'nearest-points missing in meteo configuration file') END IF !kriging-variance IF (KeyIsPresent ('kriging-variance', ini, & section = 'precipitation-daily') ) THEN krige_var = IniReadInt ('kriging-variance', ini, & section = 'precipitation-daily') CALL Catch ('info', 'PrecipitationDaily', & 'kriging variance export set to ', argument = ToString(krige_var) ) ELSE !kriging-variance missing krige_var = 0 CALL Catch ('info', 'PrecipitationDaily', & 'kriging variance export set to default (0)') END IF !kriging-anisotropy IF (KeyIsPresent ('kriging-anisotropy', ini, & section = 'precipitation-daily') ) THEN krige_anisotropy = IniReadInt ('kriging-anisotropy', ini, & section = 'precipitation-daily') CALL Catch ('info', 'PrecipitationDaily', 'kriging anisotropy set to: ', argument = ToString (krige_anisotropy) ) ELSE !kriging-anisotropy missing krige_anisotropy = 0 CALL Catch ('info', 'PrecipitationDaily', & 'kriging anisotropy set to default (0)') END IF !kriging-variogram-model IF (KeyIsPresent ('kriging-variogram-model', ini, & section = 'precipitation-daily') ) THEN krige_varmodel = IniReadInt ('kriging-variogram-model', ini, & section = 'precipitation-daily') CALL Catch ('info', 'PrecipitationDaily', 'kriging variogram model set to: ', argument = ToString (krige_varmodel) ) IF (krige_varmodel == 0 ) krige_varmodel = 5 !automatic fitting ELSE !kriging-variogram-model krige_varmodel = 2 CALL Catch ('info', 'PrecipitationDaily', & 'kriging variogram model set to default (2 exponential)') END IF !kriging-lags IF (KeyIsPresent ('kriging-lags', ini, & section = 'precipitation-daily') ) THEN krige_lags = IniReadInt ('kriging-lags', ini, & section = 'precipitation-daily') IF (krige_lags == 0) krige_lags = 15 CALL Catch ('info', 'PrecipitationDaily', 'kriging lags set to: ', argument = ToString (krige_lags) ) ELSE !kriging-variogram-model krige_lags = 15 CALL Catch ('info', 'PrecipitationDaily', & 'kriging lags set to default (15)') END IF !kriging-maxlag IF (KeyIsPresent ('kriging-maxlag', ini, & section = 'precipitation-daily') ) THEN krige_maxlag = IniReadInt ('kriging-maxlag', ini, & section = 'precipitation-daily') CALL Catch ('info', 'PrecipitationDaily', 'kriging max lag set to: ', argument = ToString (krige_maxlag) ) ELSE !kriging-maxlag missing krige_maxlag = 0 CALL Catch ('error', 'PrecipitationDaily', & 'kriging max lag set to default (0)') END IF END SELECT RETURN END SUBROUTINE SetSpecificProperties END MODULE PrecipitationDaily